Nonlinear envelope equation for broadband optical pulses in quadratic media 
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We derive a nonlinear envelope equation to describe the propagation of broadband 
optical pulses in second order nonlinear materials. The equation is first order in the 
propagation coordinate and is valid for arbitrarily wide pulse bandwidth. Our approach 
goes beyond the usual coupled wave description of x^ phenomena and provides an 
accurate modelling of the evolution of ultra-broadband pulses also when the separation 
into different coupled frequency components is not possible or not profitable. 
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The analysis of optical pulse propagation typically in- 
volves the definition of a complex envelope whose vari- 
ation is supposed to be "slow" with respect to the os- 
cillation of a carrier frequency ( "slowly varying envelope 
approximation", SVEA (l|). In the frequency domain 
this assumption is equivalent to require that the band- 
width of the envelope is narrow with respect to the car- 
rier frequency. Different works showed that it is possible 
to extend the validity of a proper generalization of the 
envelope equation (for example the "Nonlinear Envelope 
Equatio" (NEE) of Brabec and Krausz) to pulse dura- 
tion down to the single optical oscillation cycle [l-Isj and 
to the generation of high order harmonics Q. When 
second order nonlinearities are considered, the usual ap- 
proach is to write coupled equations for the separated 
frequency bands relevant for the process [si @|. How- 
ever when ultra-broadband x^ phenomena take place, 
the different frequency bands might merge, generating 
a single broad spectrum, as observed in recent experi- 
ments [lij. Obviously in these cases the coupled NEE 
description of the propagation fails due to the overlap- 
ping between different frequency bands. 

The scope of this Letter is to provide a single wave 
envelope equation to describe ultra-broadband x^ in- 
teractions. To date, such a model is not available and 
the only way to numerically describe phenomena as those 
reported in Ref. [10] is to solve directly Maxwell equa- 
tions in time domain, with an immense computational 
burden. Our equation, besides providing a powerful tool 
for analytical treatment due to its simplicity, can be eas- 
ily solved with a modest computational effort and can be 
easily generalized to include other kind of nonlinearities 
such as Kerr or Raman. 

As far the linear dispersive terms are concerned, our 
derivation of the envelope equation builds upon the 
work of Brabec and Krausz [2], that carry to a sim- 
ple model that was shown (theoretically and experimen- 
tally) to be accurate in most situations. Starting from 
Maxwell equations (written in MKS units), neglecting 
transverse dimensions (i.e considering the propagation of 



plane waves), we can obtain the 1+1D wave equation for 
the electric field E(z,t): 
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that can be written in frequency domain, by defining the 
Fourier transform F[E\{u) = E(u) = f+™ E(t)e~ iuJt dt: 
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where c is the vacuum velocity of light, eo is the vacuum 
dielectric permittivity, s(uj) = 1 + x(oo) and x(uj) is the 
linear electric susceptibility. 

We consider now the electric field E and the nonlinear 
polarization Pnl as the product of a complex envelope 



and a carrier wave: E(z,t) = A(z,t)/2e 
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P N L(z,t) = A p (z,t)/2e iuJot - if3oZ + c.c [in frequency do- 
main reads: E(z,uj) = A(z,uj-uJo)/2e~ if3oZ + A*(z, —u — 
uj )/2e if3oz }, where ujq is a reference frequency, /3q = 
Re[k(uo)] and k(u) = (u / c)^/e(uo) is the propagation 
constant. 

Particular care must be devoted to the definition of the 
complex envelope, since we do not want to put any limi- 
tation to the frequency extent of the signals. This aspect 
is commonly overlooked in literature, and it is taken for 
granted that the band of the envelope is "narrow" in some 
sense. We shall see later that for quadratically nonlin- 
ear media, a proper definition of the envelope is crucial. 
As usual in the theory of modulation [ll|, we define the 
analytic representation of the electric field: 



E(z,t) =E(z,t)+iH[E\(z,t), 
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where 
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is the Hilbert transform of the electric field (p. v. indicates 
the Cauchy principal value of the integral). The Fourier 
transform of the analytic signal reads: 
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E(z,u) = { E(z,0) if uj = 
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that is a signal that contains only the positive frequency 
content of the electric field. Due to reality of E(z,t), its 
Fourier transform has Hermitian symmetry, so that only 
the positive (or the negative) frequencies carry informa- 
tion, and we can write: 

E(z,uj) = H(z,uj) + ±fr(z,-uj), (6) 

and eventually we can define the complex electric field 
envelope as: 

A(z,t) = E(z,t)e- iuot ^ oZ , (7) 

i.e. the inverse Fourier transform of the positive fre- 
quency content of E shifted towards the low frequency 
part of the spectrum by an amount ujo- It is worth not- 
ing that no approximations on the frequency extent of 
the envelope has been done, and so supp{A(z, uj)} = 
(-cjo,+oo). 

The substitution of expressions of E(z,uu) and 
Pnl(z,oj) in Eq. ((2]), Taylor-expansion of k(uS) about 
ujo, application the slowly evolving wave approximation 
(SEWA, that is the neglect of second space derivative in 
the coordinate system moving at the group velocity at 
the reference frequency), followed by an inverse Fourier 
transform yields [l|-[3[ : 
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where D = £~ =2 ^k m {-i& ) m , k m = £( Wo ), z> = z 
and r = t — k\z is the coordinate system moving at the 
reference group velocity. 

It is worth noting that when the requirement 
| ^ ~^ ofel | << 1 is accomplished, SEWA does not ex- 
plicitly impose a limitation on pulse duration and band- 
width. Far from resonances, this requirement is fulfilled 
in the majority of parametric processes in which all waves 
propagate in the same direction. 

We now consider an instantaneous second order x^ 
nonlinearity, giving rise to the following nonlinear polar- 
ization: 



P N L(z,t)=e X i2) E(z,t) 
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It is worth noting that, due to the definition of A, the first 
(second) term in the square brackets contains only posi- 
tive (negative) frequencies, whereas the third has both. It 



is now apparent that it is impossible to separate the non- 
linear polarization in two distinct and "narrow" bands for 
the positive and negative frequencies, as common in cu- 
bic media. Moreover the neglect of the third term leads 
to totally wrong results (this term is responsible for dif- 
ference frequency generation). 

By going through the steps ((3])- (0 we can instead cor- 
rectly define the nonlinear polarization envelope: 



Ap(z,t) = P NL (z,t)e 
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Before inserting Eq. ([TO]) into Eq. (j8j), the term \A\ 2 in 
Eqs. ([9|) and ([TO]) deserves further comments, since it is 
centered around zero in frequency domain. In particu- 
lar to obtain the nonlinear polarization envelope in Eq. 
([TO]) we had to filter out the negative frequency compo- 
nents of Pjvl(^)> as done for E(uo). We not however that 
(i) A(z,uu — ujq) does not contain negative frequency by 
definition, (ii) P^l is a small perturbation to linear po- 
larization and (iii) negative frequencies cannot be phase- 
matched. It follows that the task of filtering the negative 
frequency components of \A\ 2 can be left to the prop- 
agation equation instead of having it explicitly in the 
definition of A p (z, t). In other words, when inserting Eq. 
(HDD into Eq. ®, we can write: \A\ 2 + iU[\A 2 \) « 2\A\ 2 . 
We have checked numerically the good accuracy of this 
approximation. Even if this approximation in not nec- 
essary in the numerical solution (it is straightforward to 
calculate the exact nonlinear polarization envelope in fre- 
quency domain), it is suitable to obtain a simple and 
manageable model for further analytical investigations. 
The NEE for A = A(z',t) becomes 
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or, performing derivatives, 
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Equation ([TT]) or ([T2]) constitutes the main result of this 
Letter. This nonlinear envelope equation first order in 
propagation coordinate provides a powerful means of de- 
scribing light pulse propagation in dispersive quadrati- 
cally nonlinear media. 

Starting from Eq. JTTJ it is straightforward to show 
that our equation conserves the total energy of the field, 
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i.e. ^7 \A{z' ,r)\ 2 dr = 0. It can also be shown 
that the total energy is conserved even if the non ap- 
proximated nonlinear polarization envelope A p [Eq. ([T0]) ] 
is used. 

We solved Eq. ([TT]) by split-step Fourier method ex- 
ploiting fourth order Runge-Kutta scheme for the non- 
linear step. 
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FIG. 1: Propagation of a femtosecond pulse in a PPLT crys- 
tal, a) Evolution of the field amplitude \A\ from numerical 
solution of Eq. b) Electric field amplitude at the crystal 

output, c) Comparison between coupled wave solution (dot- 
ted curves) and \A\ filtered around fundamental and second 
harmonic (solid curves). The initial pulse has gaussian shape 
and the parameters are T — 60/s, / — lOGW/cm 2 , Xi n — 
1400nra, A 2nc/oj 700nra, d 3 3 Xlt/ 2 ^0.6pm/V 

In order to show the validity of our equation, we 
simulated the propagation of a femtosecond pulse in a 
L = bmm long periodically poled lithium tantalate sam- 
ple (PPLT). To model the refractive index dispersion 
we employed a Sellmeier model fitted from experimen- 

I li (2) 

tal data [12j and nonlinear coefficient is dss = Xlt/ 2 = 
10. 6pm/V. In the numerical code we inserted the ex- 
act dispersion relation k(uj). We assumed a first or- 
der quasi phase matching (QPM) grating, with a period 
A = llAjim. We thus allowed a periodic variation of the 
nonlinear coefflcent X^ = X <y2 \ z ) — 2/7rx^e* 27r / A2; +c.c.. 
We injected a T = 60fs FWHM long gaussian pulse, 
centered around 1400nm, with / = lOGW/cm 2 peak in- 
tensity. The corresponding residual phase mismatch is 
Ak = 2k(u in ) - k(2u in ) = lOOOOra -1 , where UJi n IS the 
carrier frequency of the input pulse. In the simulation 



we set the reference frequency ujq to be equal to the sec- 
ond harmonic of the input pulse: in this way the second 
harmonic is stationary in the reference frame (z',t). 
Figure QJi) shows the evolution of the electric field enve- 
lope amplitude \A\ from numerical solution of Eq. ([TT]) . 
We can see the typical scenario of the propagation of 
femtosecond pulses in highly group velocity mismatched 
(GVM) process: the fundamental frequency (FF) pulse 
generates its second harmonic (SH) during propagation, 
and the generated SH pulse has the typical shape of a ini- 
tial peak followed by a long tail, whose duration is fixed 
by the product between GVM and crystal length. Figure 
[Ud) shows the electric field envelope amplitude at the end 
of the crystal. It can be seen a peak that corresponds to 
the faster frequency components located around FF, fol- 
lowed by a long tail that ends with a second lower peak. 
This long pulse corresponds to the generated SH compo- 
nents. The SH pulse is smooth, indicating that no beat- 
ing with eventual FF components is present. Whereas 
in the residual FF pulses centered around r ~ — 1.6ps, 
there is a clear fast oscillation, indicating that FF and 
SH components are superimposed. 

To test the results, we simulated the same set-up with a 
standard coupled wave model [13], by inserting the val- 
ues of first and second order dispersion evaluated at FF 
and SH. To compare the results we filtered A around FF 
and SH. Figure Ht) shows the electric field amplitudes 
at the end of the crystal. The results of the two mod- 
els are practically indistinguishable. It is worth noting 
that this simulation shows the validity of Eq. ([TT]) over 
a bandwidth of uoq. 

As a second example we consider the propagation of 
a femtosecond pulse into a highly mismatched period- 
ically poled lithium niobate (PPLN) sample, that was 
demonstrated experimentally to generate an octave span- 
ning supercontinuum spectral broadening [To| . To model 
the refractive index dispersion we employed a Sellmeier 
model fitted from experimental data [14[ and nonlinear 

(2) 

coefficient is ^33 = Xln/ 2 ~ 27pra/V. In the numerical 
code we inserted the exact dispersion relation k(uo). We 
assumed a QPM grating with a period A = 30/ira (phase 
matched for second harmonic generation at around 2\im 
fundamental wavelength). We included higher order 
QPM terms, since the huge bandwidth can phase match 
different spatial harmonics. We injected a T = 50 fs 
FWHM long gaussian pulse, centered around 1580nra, 
with / = lbGW/cm 2 peak intensity. In the simulation 
we set the reference wavelength Ao = 700nm. Figure [2^i) 
shows the evolution of the spectrum during the propaga- 
tion into a L = 7mm crystal. We can see a consistent 
broadening and redshift of the FF part of the spectrum 
that, at the end of the crystal, reaches an octave-spanning 
bandwidth from 1200nm to 3000nm. We can also see 
the generation of spectral components at the second and 
third harmonics. At the second harmonic the spectrum 
initially broadens and has an evolution ruled by highly 
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FIG. 2: Propagation of a femtosecond pulse into a PPLN 
sample, a) Evolution of the power spectrum (in dB) from 
numerical solution of Eq. ([TT]) . b) Power specrum at the 
crystal output in the visible and NIR range, c) Power specrum 
at the crystal output in the infrared range. The division into 
separate spectral region is made to facilitate the comparison 
with experimental data [10]. The initial pulse has gaussian 
shape and the parameters are T = 50/s, / = 15GW/cra 2 , 
Xin — 1580nra, Ao = 2ttc/ujo — 700nra, 0I33 — = 
27pm /V 

mismatched SHG. When the FF broadening reaches the 
first order quasi phase matching wavelength at around 
2/xra, the more efficient conversion process generates a 
spike at around 1/ira. Figure [2] b) shows the visible and 
the near infrared (NIR) part of the spectrum at the crys- 
tal output. We can see a broadband second and third 
harmonic of the broadened laser spectrum, and the pres- 
ence of some spikes given by the quasi phase matching 
of high order spatial harmonics of the grating. We ver- 
ified that the two spikes at the third harmonic corre- 
spond to the third and fifth order QPM for the process 
u + 2u — )> 3cj. We can also see a spectral overlap between 
the harmonics of the broadened laser spectrum, that can 
be exploited to achieve carrier-envelope-offset phase slip 
stabilization [10j, that is of paramount importance for 
frequency metrology applications. 

Figure [2t) shows the infrared spectrum at the output. 
This spectrum exhibits more than an octave spanning be- 
tween 1300nra and 3000nra at the —40dB spectral power 
lever with respect to the peak power level. The spec- 
tral components near the zero GVM wavelength around 



3000nm are generated more efficiently. 
All the features described above compares surprisi ngly 
well with the experimental results of Langrock et al. [10j , 
even if we simulate a slightly different environment. In 
fact we use a bulk PPLN sample and not a RPE PPLN 
waveguide. The effect of waveguide is to slightly modify 
the power levels and the crystal dispersion: a detailed 
simulation of the real set-up is out of the scope of this 
Letter. It is worth noting that numerical modelling of 
such phenomena without our model is an irksome job 
since (i) time domain Maxwell equation solvers require 
a prohibitive computational effort and (ii) coupled wave 
approaches cannot be used in the presence of overlapping 
among frequency bands of different field components. 

In conclusion we have derived a robust nonlinear en- 
velope equation describing the propagation in dispersive 
quadratic materials. Thanks to a proper formal definition 
of the complex envelope, it is possible to treat pulses of 
arbitrary frequency content. A proper definition of enve- 
lope is crucial for second order nonlinearities, due to the 
generation of frequency components around zero. Com- 
putationally it is possible to accurately evolve optical 
pulses of arbitrarily wide band over a meter scale phys- 
ical distance, which is a few order of magnitude longer 
than those accessible by Maxwell equation solvers. 
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